home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 21 / Cream of the Crop 21 (Terry Blount) (October 1996).iso / os2 / e33el2.zip / emacs / 19.33 / lisp / float.el < prev    next >
Lisp/Scheme  |  1996-01-20  |  16KB  |  459 lines
  1. ;;; float.el --- floating point arithmetic package.
  2.  
  3. ;; Copyright (C) 1986 Free Software Foundation, Inc.
  4.  
  5. ;; Author: Bill Rosenblatt
  6. ;; Maintainer: FSF
  7. ;; Keywords: extensions
  8.  
  9. ;; This file is part of GNU Emacs.
  10.  
  11. ;; GNU Emacs is free software; you can redistribute it and/or modify
  12. ;; it under the terms of the GNU General Public License as published by
  13. ;; the Free Software Foundation; either version 2, or (at your option)
  14. ;; any later version.
  15.  
  16. ;; GNU Emacs is distributed in the hope that it will be useful,
  17. ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19. ;; GNU General Public License for more details.
  20.  
  21. ;; You should have received a copy of the GNU General Public License
  22. ;; along with GNU Emacs; see the file COPYING.  If not, write to the
  23. ;; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  24. ;; Boston, MA 02111-1307, USA.
  25.  
  26. ;;; Commentary:
  27.  
  28. ;; Floating point numbers are represented by dot-pairs (mant . exp)
  29. ;; where mant is the 24-bit signed integral mantissa and exp is the
  30. ;; base 2 exponent.
  31. ;;
  32. ;; Emacs LISP supports a 24-bit signed integer data type, which has a
  33. ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
  34. ;; This gives six significant decimal digit accuracy.  Exponents can
  35. ;; be anything in the range -(2**23) to +(2**23)-1.
  36. ;;
  37. ;; User interface:
  38. ;; function f converts from integer to floating point
  39. ;; function string-to-float converts from string to floating point
  40. ;; function fint converts a floating point to integer (with truncation)
  41. ;; function float-to-string converts from floating point to string
  42. ;;                   
  43. ;; Caveats:
  44. ;; -  Exponents outside of the range of +/-100 or so will cause certain 
  45. ;;    functions (especially conversion routines) to take forever.
  46. ;; -  Very little checking is done for fixed point overflow/underflow.
  47. ;; -  No checking is done for over/underflow of the exponent
  48. ;;    (hardly necessary when exponent can be 2**23).
  49. ;; 
  50. ;;
  51. ;; Bill Rosenblatt
  52. ;; June 20, 1986
  53. ;;
  54.  
  55. ;;; Code:
  56.  
  57. ;; fundamental implementation constants
  58. (defconst exp-base 2
  59.   "Base of exponent in this floating point representation.")
  60.  
  61. (defconst mantissa-bits 24
  62.   "Number of significant bits in this floating point representation.")
  63.  
  64. (defconst decimal-digits 6
  65.   "Number of decimal digits expected to be accurate.")
  66.  
  67. (defconst expt-digits 2
  68.   "Maximum permitted digits in a scientific notation exponent.")
  69.  
  70. ;; other constants
  71. (defconst maxbit (1- mantissa-bits)
  72.   "Number of highest bit")
  73.  
  74. (defconst mantissa-maxval (1- (ash 1 maxbit))
  75.   "Maximum permissible value of mantissa")
  76.  
  77. (defconst mantissa-minval (ash 1 maxbit)
  78.   "Minimum permissible value of mantissa")
  79.  
  80. (defconst floating-point-regexp
  81.   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
  82. \\(\\.\\([0-9]*\\)\\|\\)\
  83. \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
  84.   "Regular expression to match floating point numbers.  Extract matches:
  85. 1 - minus sign
  86. 2 - integer part
  87. 4 - fractional part
  88. 8 - minus sign for power of ten
  89. 9 - power of ten
  90. ")
  91.  
  92. (defconst high-bit-mask (ash 1 maxbit)
  93.   "Masks all bits except the high-order (sign) bit.")
  94.  
  95. (defconst second-bit-mask (ash 1 (1- maxbit))
  96.   "Masks all bits except the highest-order magnitude bit")
  97.  
  98. ;; various useful floating point constants
  99. (setq _f0 '(0 . 1))
  100.  
  101. (setq _f1/2 '(4194304 . -23))
  102.  
  103. (setq _f1 '(4194304 . -22))
  104.  
  105. (setq _f10 '(5242880 . -19))
  106.  
  107. ;; support for decimal conversion routines
  108. (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
  109. (aset powers-of-10 1 _f10)
  110. (aset powers-of-10 2 '(6553600 . -16))
  111. (aset powers-of-10 3 '(8192000 . -13))
  112. (aset powers-of-10 4 '(5120000 . -9))
  113. (aset powers-of-10 5 '(6400000 . -6))
  114. (aset powers-of-10 6 '(8000000 . -3))
  115.  
  116. (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
  117.       highest-power-of-10 (aref powers-of-10 decimal-digits))
  118.  
  119. (defun fashl (fnum)            ; floating-point arithmetic shift left
  120.   (cons (ash (car fnum) 1) (1- (cdr fnum))))
  121.  
  122. (defun fashr (fnum)            ; floating point arithmetic shift right
  123.   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
  124.  
  125. (defun normalize (fnum)
  126.   (if (> (car fnum) 0)            ; make sure next-to-highest bit is set
  127.       (while (zerop (logand (car fnum) second-bit-mask))
  128.     (setq fnum (fashl fnum)))
  129.     (if (< (car fnum) 0)        ; make sure highest bit is set
  130.     (while (zerop (logand (car fnum) high-bit-mask))
  131.       (setq fnum (fashl fnum)))
  132.       (setq fnum _f0)))            ; "standard 0"
  133.   fnum)
  134.       
  135. (defun abs (n)                ; integer absolute value
  136.   (if (>= n 0) n (- n)))
  137.  
  138. (defun fabs (fnum)            ; re-normalize after taking abs value
  139.   (normalize (cons (abs (car fnum)) (cdr fnum))))
  140.  
  141. (defun xor (a b)            ; logical exclusive or
  142.   (and (or a b) (not (and a b))))
  143.  
  144. (defun same-sign (a b)            ; two f-p numbers have same sign?
  145.   (not (xor (natnump (car a)) (natnump (car b)))))
  146.  
  147. (defun extract-match (str i)        ; used after string-match
  148.   (condition-case ()
  149.       (substring str (match-beginning i) (match-end i))
  150.     (error "")))
  151.  
  152. ;; support for the multiplication function
  153. (setq halfword-bits (/ mantissa-bits 2)    ; bits in a halfword
  154.       masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
  155.       maskhi (lognot masklo)        ; isolate the upper halfword
  156.       round-limit (ash 1 (/ halfword-bits 2)))
  157.  
  158. (defun hihalf (n)            ; return high halfword, shifted down
  159.   (ash (logand n maskhi) (- halfword-bits)))
  160.  
  161. (defun lohalf (n)            ; return low halfword
  162.   (logand n masklo))
  163.  
  164. ;; Visible functions
  165.  
  166. ;; Arithmetic functions
  167. (defun f+ (a1 a2)
  168.   "Returns the sum of two floating point numbers."
  169.   (let ((f1 (fmax a1 a2))
  170.     (f2 (fmin a1 a2)))
  171.     (if (same-sign a1 a2)
  172.     (setq f1 (fashr f1)        ; shift right to avoid overflow
  173.           f2 (fashr f2)))
  174.     (normalize
  175.      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
  176.        (cdr f1)))))
  177.  
  178. (defun f- (a1 &optional a2)        ; unary or binary minus
  179.   "Returns the difference of two floating point numbers."
  180.   (if a2
  181.       (f+ a1 (f- a2))
  182.     (normalize (cons (- (car a1)) (cdr a1)))))
  183.  
  184. (defun f* (a1 a2)            ; multiply in halfword chunks
  185.   "Returns the product of two floating point numbers."
  186.   (let* ((i1 (car (fabs a1)))
  187.      (i2 (car (fabs a2)))
  188.      (sign (not (same-sign a1 a2)))
  189.      (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
  190.             (lohalf (* (hihalf i1) (lohalf i2)))
  191.             (lohalf (* (lohalf i1) (hihalf i2)))))
  192.      (prodhi (+ (* (hihalf i1) (hihalf i2))
  193.             (hihalf (* (hihalf i1) (lohalf i2)))
  194.             (hihalf (* (lohalf i1) (hihalf i2)))
  195.             (hihalf prodlo))))
  196.     (if (> (lohalf prodlo) round-limit)
  197.     (setq prodhi (1+ prodhi)))    ; round off truncated bits
  198.     (normalize
  199.      (cons (if sign (- prodhi) prodhi)
  200.        (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
  201.  
  202. (defun f/ (a1 a2)            ; SLOW subtract-and-shift algorithm
  203.   "Returns the quotient of two floating point numbers."
  204.   (if (zerop (car a2))            ; if divide by 0
  205.       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
  206.     (let ((bits (1- maxbit))
  207.       (quotient 0) 
  208.       (dividend (car (fabs a1)))
  209.       (divisor (car (fabs a2)))
  210.       (sign (not (same-sign a1 a2))))
  211.       (while (natnump bits)
  212.     (if (< (- dividend divisor) 0)
  213.         (setq quotient (ash quotient 1))
  214.       (setq quotient (1+ (ash quotient 1))
  215.         dividend (- dividend divisor)))
  216.     (setq dividend (ash dividend 1)
  217.           bits (1- bits)))
  218.       (normalize
  219.        (cons (if sign (- quotient) quotient)
  220.          (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
  221.   
  222. (defun f% (a1 a2)
  223.   "Returns the remainder of first floating point number divided by second."
  224.   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
  225.       
  226.  
  227. ;; Comparison functions
  228. (defun f= (a1 a2)
  229.   "Returns t if two floating point numbers are equal, nil otherwise."
  230.   (equal a1 a2))
  231.  
  232. (defun f> (a1 a2)
  233.   "Returns t if first floating point number is greater than second,
  234. nil otherwise."
  235.   (cond ((and (natnump (car a1)) (< (car a2) 0)) 
  236.      t)                ; a1 nonnegative, a2 negative
  237.     ((and (> (car a1) 0) (<= (car a2) 0))
  238.      t)                ; a1 positive, a2 nonpositive
  239.     ((and (<= (car a1) 0) (natnump (car a2)))
  240.      nil)                ; a1 nonpos, a2 non